package scu.maqiang.fes;

import java.util.Arrays;
import java.util.HashSet;
import java.util.Map;

import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

import static scu.maqiang.numeric.Constants.*;

/**
 * 二维线性三角形有限元空间, 对于二维情况下热传导, 电磁学等问题中整体矩阵进行组装
 */
public class FES2T31 extends FES {
    public FES2T31(Mesh2T3 mesh2D) {
        super(mesh2D, 1);
        eleStiff = FES2T31::elementStiff;
        eleMass = FES2T31::elementMass;
        eleSource = FES2T31::elementSource;
        eleFlux = FES2T31::elementFlux;
    }

    public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
        return getValue0(u, xyz, eleIdx, eleVal, CG::pointInT3);
    }

    public double getValue(double[] u, double[] xyz) {
        return getValue0(u, xyz, CG::pointInT3);
    }

    public void getValue(double[] u, double[] xyz, DoubleIntegerPair data, int[] eleIdx, double[] eleVal) {
        getValue0(u, xyz, data, eleIdx, eleVal, CG::pointInT3);
    }

    public SRMatrix assembleStiff(ScalarFunc[] coefFunc, double[] param, BVPType tp) {
        SRMatrix A = new SRMatrix(ndof);
        assembleGlobalMatrix(coefFunc, param, FES2T31::elementStiff, tp, A);
        return A;
    }

    public void assembleStiff(double[] constCoef, BVPType tp, SCMatrix A, int... label) {
        assembleGlobalMatrix(constCoef, FES2T31::elementStiff, tp, A, label);
    }

//    public void assembleNLStiff(double[] constCoef, double[][] vecCoef, BVPType tp, SRMatrix A) {
//        assembleGlobalMatrix(constCoef, vecCoef, FES2T31::elementStiff, tp, A);
//    }


    /**
     * 计算单元热刚度矩阵
     *
     * @param xy        单元各结点坐标
     * @param kxx       热传导系数
     * @param tp        矩阵类型
     * @param eleMatrix 存储单元刚度矩阵
     */
    public static void elementStiff(int[] ele, double[][] xy, double[][] kxx, BVPType tp, double[][] eleMatrix) {

        double[][] grad = new double[3][2];
        double area = CG.gradient2DT3(xy, grad);
        double[] dx = MVO.col(grad, 0);
        double[] dy = MVO.col(grad, 1);
        double dxU, dyU;
        //System.out.println(area);
        double[] x = MVO.col(xy, 0);
        double Af;
        double c11, c22;
        double[][] temp;
        double[] N = new double[3];
        double[][] gradN = new double[3][2];
        double[] coef0_11, coefX_11, coefX2_11, coefX3_11, coefU_11, coefU2_11, coefU3_11;
        double[] coef0_22, coefX_22, coefX2_22, coefX3_22, coefU_22, coefU2_22, coefU3_22;
        double[] coefdX_11, coefdX2_11, coefdX3_11, coefdY_11, coefdY2_11, coefdY3_11;
        double[] coefdX_22, coefdX2_22, coefdX3_22, coefdY_22, coefdY2_22, coefdY3_22;
        double[] u;
        double weight;
        MVO.fill(eleMatrix, 0.0);
        switch (tp) {
            case COMMON:
                Af = area * MVO.average(kxx[0]);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1]) * Af;
                    }
                }
                break;
            case AXISYMMETRIC:
                Af = area * (MVO.sum(kxx[0]) * MVO.sum(x) + MVO.dot_product(kxx[0], x)) / 12.0;
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1]) * Af;
                    }
                }
                break;
            case CONSTITUTE_MATRIX_COMMON:
                temp = MVO.matmul(grad, kxx);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = 0.0;
                        for (int k = 0; k < 2; k++) {
                            eleMatrix[i][j] += temp[i][k] * grad[j][k] * area;
                        }
                    }
                }

                break;
            case CONSTITUTE_MATRIX_AXIS:
                temp = MVO.matmul(grad, kxx);
                Af = MVO.average(x);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = 0.0;
                        for (int k = 0; k < 2; k++) {
                            eleMatrix[i][j] += temp[i][k] * grad[j][k] * Af;
                        }
                    }
                }
                break;
            case NONLINEAR_POLYNORMIAL_2:
                //a11 + f11 u + g11 u^2 + h11 u^3
                //a22 + f22 u + g22 u^2 + h22 u^3
                coef0_11 = kxx[0];
                coefU_11 = kxx[1];
                coefU2_11 = kxx[2];
                coefU3_11 = kxx[3];
                coef0_22 = kxx[4];
                coefU_22 = kxx[5];
                coefU2_22 = kxx[6];
                coefU3_22 = kxx[7];
                u = kxx[8];
                dxU = MVO.dot_product(u, dx);
                dyU = MVO.dot_product(u, dy);
                for (int i = 0; i < 7; i++) {
                    c11 = 0;
                    c22 = 0;
                    N = Constants.TP5[i].clone();
                    double meanU = MVO.dot_product(u, N);
                    c11 += MVO.dot_product(coefU_11, N);
                    c11 += 2 * meanU * MVO.dot_product(coefU2_11, N);
                    c11 += 3 * meanU * meanU * MVO.dot_product(coefU3_11, N);
                    c22 += MVO.dot_product(coefU_22, N);
                    c22 += 2 * meanU * MVO.dot_product(coefU2_22, N);
                    c22 += 3 * meanU * meanU * MVO.dot_product(coefU3_22, N);
                    weight = area * Constants.TW5[i];
                    for (int si = 0; si < 3; si++) {
                        for (int sj = 0; sj < 3; sj++) {
                            eleMatrix[si][sj] += (c11 * dxU * grad[si][0] * N[sj] + c22 * dyU * grad[si][1] * N[sj]) * weight;
                        }
                    }
                }
            case NONLINEAR_POLYNORMIAL_1:
                //a11 + f11 u + g11 u^2 + h11 u^3
                //a22 + f22 u + g22 u^2 + h22 u^3
                coef0_11 = kxx[0];
                coefU_11 = kxx[1];
                coefU2_11 = kxx[2];
                coefU3_11 = kxx[3];
                coef0_22 = kxx[4];
                coefU_22 = kxx[5];
                coefU2_22 = kxx[6];
                coefU3_22 = kxx[7];
                u = kxx[8];
                for (int i = 0; i < 3; i++) {
                    c11 = 0;
                    c22 = 0;
                    N = Constants.TP2[i].clone();
                    c11 += MVO.dot_product(coef0_11, N);
                    c22 += MVO.dot_product(coef0_22, N);
                    double meanU = MVO.dot_product(u, N);
                    c11 += MVO.dot_product(coefU_11, N) * meanU;
                    c22 += MVO.dot_product(coefU_22, N) * meanU;
                    c11 += MVO.dot_product(coefU2_11, N) * meanU * meanU;
                    c22 += MVO.dot_product(coefU2_22, N) * meanU * meanU;
                    c11 += MVO.dot_product(coefU3_11, N) * meanU * meanU * meanU;
                    c22 += MVO.dot_product(coefU3_22, N) * meanU * meanU * meanU;
                    //System.out.println("c11 = " + c11 + "\tc22 = " + c22);
                    weight = area * Constants.TW2[i];
                    for (int si = 0; si < 3; si++) {
                        for (int sj = 0; sj < 3; sj++) {
                            eleMatrix[si][sj] += (c11 * grad[si][0] * grad[sj][0] + c22 * grad[si][1] * grad[sj][1]) * weight;
                        }
                    }
                }
                break;
            case NONLINEAR_DERIVATIVE_2:
                //a11 + b11 du/dx + c11 (du/dx)^2 + d11 (du/dx)^3 + f11 du/dy + g11 (du/dy)^2 + h11 (du/dy)^3
                //a22 + b22 du/dx + d22 (du/dx)^2 + d22 (du/dx)^3 + f22 du/dy + g22 (du/dy)^2 + h22 (du/dy)^3
                double c12, c21;
                coef0_11 = kxx[0];
                coefdX_11 = kxx[1];
                coefdX2_11 = kxx[2];
                coefdX3_11 = kxx[3];
                coefdY_11 = kxx[4];
                coefdY2_11 = kxx[5];
                coefdY3_11 = kxx[6];
                coef0_22 = kxx[7];
                coefdX_22 = kxx[8];
                coefdX2_22 = kxx[9];
                coefdX3_22 = kxx[10];
                coefdY_22 = kxx[11];
                coefdY2_22 = kxx[12];
                coefdY3_22 = kxx[13];
                u = kxx[14];
                dxU = MVO.dot_product(dx, u);
                dyU = MVO.dot_product(dy, u);
                c11 = (MVO.average(coefdX_11) + (2 * MVO.average(coefdX2_11) + 3 * MVO.average(coefdX3_11) * dxU) * dxU) * dxU;
                c12 = (MVO.average(coefdY_11) + (2 * MVO.average(coefdY2_11) + 3 * MVO.average(coefdY3_11) * dyU) * dyU) * dxU;
                c21 = (MVO.average(coefdX_22) + (2 * MVO.average(coefdX2_22) + 3 * MVO.average(coefdX3_22) * dxU) * dxU) * dyU;
                c22 = (MVO.average(coefdY_22) + (2 * MVO.average(coefdY2_22) + 3 * MVO.average(coefdY3_22) * dyU) * dyU) * dyU;
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = (c11 * grad[i][0] * grad[j][0] + c12 * grad[i][0] * grad[j][1]
                                          +c21 * grad[i][1] * grad[j][0] + c22 * grad[i][1] * grad[j][1]) * area;
                    }
                }
            case NONLINEAR_DERIVATIVE_1:
                coef0_11   = kxx[0];
                coefdX_11  = kxx[1];
                coefdX2_11 = kxx[2];
                coefdX3_11 = kxx[3];
                coefdY_11  = kxx[4];
                coefdY2_11 = kxx[5];
                coefdY3_11 = kxx[6];
                coef0_22   = kxx[7];
                coefdX_22  = kxx[8];
                coefdX2_22 = kxx[9];
                coefdX3_22 = kxx[10];
                coefdY_22  = kxx[11];
                coefdY2_22 = kxx[12];
                coefdY3_22 = kxx[13];
                u = kxx[14];
                dxU = MVO.dot_product(u, dx);
                dyU = MVO.dot_product(u, dy);
                c11 = MVO.average(coef0_11);
                c11 += MVO.average(coefdX_11)  * dxU;
                c11 += MVO.average(coefdX2_11) * dxU * dxU;
                c11 += MVO.average(coefdX3_11) * dxU * dxU * dxU;
                c11 += MVO.average(coefdY_11)  * dyU;
                c11 += MVO.average(coefdY2_11) * dyU * dyU;
                c11 += MVO.average(coefdY3_11) * dyU * dyU * dyU;
                c22 = MVO.average(coef0_22);
                c22 += MVO.average(coefdX_22)  * dxU;
                c22 += MVO.average(coefdX2_22) * dxU * dxU;
                c22 += MVO.average(coefdX3_22) * dxU * dxU * dxU;
                c22 += MVO.average(coefdY_22)  * dyU;
                c22 += MVO.average(coefdY2_22) * dyU * dyU;
                c22 += MVO.average(coefdY3_22) * dyU * dyU * dyU;
                for (int si = 0; si < 3; si++) {
                    for (int sj = 0; sj < 3; sj++) {
                        eleMatrix[si][sj] = (c11 * grad[si][0] * grad[sj][0] + c22 * grad[si][1] * grad[sj][1]) * area;
                    }
                }
                break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
        // System.out.println(MVO.toString(kxx[0]));
        // System.out.println(Arrays.toString(kxx[0]));
        // System.out.println(MVO.toString(eleMatrix));
//
//
// 		System.out.println(MVO.toString(eleMatrix));
//		System.exit(1);
    }

    public void assembleStiff(MatrixCFunc coefFunc, double[][][] param, BVPType tp, SCMatrix A, int... label) {
        assembleGlobalMatrix(coefFunc, param, FES2T31::elementCStiff, tp, A, label);
    }

//	assembleGlobalMatrix(MatrixCFunc matrixCoefFunc, double[] param, EleMatCFunc computeEleMa, BVPType tp,
//						 SCMatrix A, int... domainLabel) {

    public static void elementCStiff(int[] ele, double[][] xy, Complex[][] kxx, BVPType tp, Complex[][] eleMatrix) {

        double[][] grad = new double[3][2];
        double area = CG.gradient2DT3(xy, grad);
        double[] x = MVO.col(xy, 0);
        Complex[][] temp;
        Complex Af;
        switch (tp) {
            case COMMON:
                Af = Complex.average(kxx[0]).mul(area);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        double innertemp = grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1];
                        eleMatrix[i][j].real = Af.real * innertemp;
                        eleMatrix[i][j].image = Af.image * innertemp;
                    }
                }
                break;
            case AXISYMMETRIC:
                Af = ((Complex.sum(kxx[0]).mul(MVO.sum(x))).add(Complex.dot_product(kxx[0], x))).mul(area / 12.0);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
//						eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1]) * Af;
                        double innertemp = grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1];
                        eleMatrix[i][j].real = Af.real * innertemp;
                        eleMatrix[i][j].image = Af.image * innertemp;
                    }
                }
                break;
            case CONSTITUTE_MATRIX_COMMON:
                temp = Complex.mul(grad, kxx);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j].real = 0.0;
                        eleMatrix[i][j].image = 0.0;
                        for (int k = 0; k < 2; k++) {
                            eleMatrix[i][j].real += temp[i][k].real * grad[j][k] * area;
                            eleMatrix[i][j].image += temp[i][k].image * grad[j][k] * area;

                        }
                        //eleMatrix[j][i].real = eleMatrix[i][j].real;
                        //eleMatrix[j][i].image = eleMatrix[i][j].image;
                    }
                }

                break;
            case CONSTITUTE_MATRIX_AXIS:
                temp = Complex.mul(grad, kxx);
                double AAf = MVO.average(x);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j].real = 0.0;
                        eleMatrix[i][j].image = 0.0;
                        for (int k = 0; k < 2; k++) {
                            eleMatrix[i][j].real += temp[i][k].real * grad[j][k] * AAf;
                            eleMatrix[i][j].image += temp[i][k].image * grad[j][k] * AAf;
                        }
                    }
                }
                break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
//		System.out.println(Complex.toString(eleMatrix[0]));
//		System.out.println(Complex.toString(eleMatrix[1]));
//		System.out.println(Complex.toString(eleMatrix[2]));
//		System.exit(1);
        // System.out.println(Arrays.toString(kxx[0]));
        // System.out.println(MVO.toString(eleMatrix));
    }

    public void assemblePMLStiffMass(double mur, double k0, double coefM, BVPType tpA, BVPType tpM, SCMatrix A, int interfaceBoundaryLabel, int outerBoundaryLabel, int... pmlDomainLabel) {
        //Pair<Integer, double[]>[] PMLnodes = mesh.getNodeCoordinates(pmlDomainLabel);
        Map<Integer, Complex[]> pmlCoordinates = mesh.setPMLCoordinates(k0, new int[]{interfaceBoundaryLabel}, new int[]{outerBoundaryLabel}, pmlDomainLabel);
        Complex[][] coord = new Complex[mesh.getnPerElement()][];
        double[][] kxx = new double[1][mesh.getnPerElement()];
        double[][] ma = new double[1][mesh.getnPerElement()];
        for (int i = 0; i < mesh.getNt(); i++) {
            int[] ele = mesh.getElement(i);
            int eleLabel = mesh.getElementLabel(i);
            for (int j = 0; j < pmlDomainLabel.length; j++) {
                if (eleLabel == pmlDomainLabel[j]) {
//					System.out.println(i + "\t" + Arrays.toString(ele));
                    for (int k = 0; k < mesh.getnPerElement(); k++) {
                        coord[k] = pmlCoordinates.get(ele[k]);
                        kxx[0][k] = 1 / mur;
                        ma[0][k] = coefM;
                    }

//					System.out.println(Arrays.toString(kxx[0]));
//					System.out.println(Complex.toString(coord[0]));
//					System.out.println(Complex.toString(coord[1]));
//					System.out.println(Complex.toString(coord[2]));
                    Complex[][] eleMatrix = Complex.zeros(mesh.getnPerElement(), mesh.getnPerElement());
                    Complex[][] eleMatrix2 = Complex.zeros(mesh.getnPerElement(), mesh.getnPerElement());
                    elementPMLStiff(ele, coord, kxx, tpA, eleMatrix);
                    elementPMLMass(ele, coord, ma, tpM, eleMatrix2);
                    A.assemble(getDof(ele), eleMatrix);
                    A.assemble(getDof(ele), eleMatrix2);
                }
            }
        }
    }

    /**
     * 计算单元热刚度矩阵
     *
     * @param xy        单元各结点坐标
     * @param kxx       热传导系数
     * @param tp        矩阵类型
     * @param eleMatrix 存储单元刚度矩阵
     */
    public static void elementPMLStiff(int[] ele, Complex[][] xy, double[][] kxx, BVPType tp, Complex[][] eleMatrix) {

        Complex[][] grad = new Complex[3][2];
        Complex area = CG.gradient2DT3(xy, grad);
        Complex[] x = Complex.Col(xy, 0);
        Complex Af;
        switch (tp) {
            case COMMON:
                Af = area.mul(MVO.average(kxx[0]));
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = ((grad[i][0].mul(grad[j][0])).add(grad[i][1].mul(grad[j][1]))).mul(Af);
                    }
                }
                break;
//			case AXISYMMETRIC:
//				Af = area * (MVO.sum(kxx[0]) * MVO.sum(x) + MVO.dot_product(kxx[0], x)) / 12.0;
//				for (int i = 0; i < 3; i++) {
//					for (int j = 0; j < 3; j++) {
//						eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1]) * Af;
//					}
//				}
//				break;
            case CONSTITUTE_MATRIX_COMMON:
                Complex[][] temp = Complex.mul(grad, kxx);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j].real = 0.0;
                        eleMatrix[i][j].image = 0.0;
                        for (int k = 0; k < 2; k++) {
                            Complex temp2 = (temp[i][k].mul(grad[j][k])).mul(area);
                            eleMatrix[i][j].real += temp2.real;
                            eleMatrix[i][j].image += temp2.image;
                        }
                    }
                }

                break;
//			case CONSTITUTE_MATRIX_AXIS:
//				temp = MVO.matmul(grad, kxx);
//				Af = MVO.average(x);
//				for (int i = 0; i < 3; i++) {
//					for (int j = 0; j < 3; j++) {
//						eleMatrix[i][j] = 0.0;
//						for (int k = 0; k < 2; k++) {
//							eleMatrix[i][j] += temp[i][k] * grad[j][k] * Af;
//						}
//					}
//				}
//				break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
        // System.out.println(MVO.toString(kxx[0]));
        // System.out.println(Arrays.toString(kxx[0]));
        // System.out.println(MVO.toString(eleMatrix));
    }

    /**
     * 计算单元热刚度矩阵
     *
     * @param xy        单元各结点坐标
     * @param kxx       热传导系数
     * @param tp        矩阵类型
     * @param eleMatrix 存储单元刚度矩阵
     */
    public static void elementPMLStiff(Complex[][] xy, Complex[][] kxx, BVPType tp, Complex[][] eleMatrix) {

        Complex[][] grad = new Complex[3][2];
        Complex area = CG.gradient2DT3(xy, grad);
        Complex[] x = Complex.Col(xy, 0);
        Complex Af;
        switch (tp) {
            case COMMON:
                Af = area.mul(Complex.average(kxx[0]));
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j] = ((grad[i][0].mul(grad[j][0])).add(grad[i][1].mul(grad[j][1]))).mul(Af);
                    }
                }
                break;
//			case AXISYMMETRIC:
//				Af = area * (MVO.sum(kxx[0]) * MVO.sum(x) + MVO.dot_product(kxx[0], x)) / 12.0;
//				for (int i = 0; i < 3; i++) {
//					for (int j = 0; j < 3; j++) {
//						eleMatrix[i][j] = (grad[i][0] * grad[j][0] + grad[i][1] * grad[j][1]) * Af;
//					}
//				}
//				break;
            case CONSTITUTE_MATRIX_COMMON:
                Complex[][] temp = Complex.mul(grad, kxx);
                for (int i = 0; i < 3; i++) {
                    for (int j = 0; j < 3; j++) {
                        eleMatrix[i][j].real = 0.0;
                        eleMatrix[i][j].image = 0.0;
                        for (int k = 0; k < 2; k++) {
                            Complex temp2 = (temp[i][k].mul(grad[j][k])).mul(area);
                            eleMatrix[i][j].real += temp2.real;
                            eleMatrix[i][j].image += temp2.image;
                        }
                    }
                }

                break;
//			case CONSTITUTE_MATRIX_AXIS:
//				temp = MVO.matmul(grad, kxx);
//				Af = MVO.average(x);
//				for (int i = 0; i < 3; i++) {
//					for (int j = 0; j < 3; j++) {
//						eleMatrix[i][j] = 0.0;
//						for (int k = 0; k < 2; k++) {
//							eleMatrix[i][j] += temp[i][k] * grad[j][k] * Af;
//						}
//					}
//				}
//				break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
        // System.out.println(MVO.toString(kxx[0]));
        // System.out.println(Arrays.toString(kxx[0]));
        // System.out.println(MVO.toString(eleMatrix));
    }

    public void assembleMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A) {
        assembleGlobalMatrix(coefFunc, param, FES2T31::elementMass, tp, A);
    }

    public void assembleMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... label) {
        assembleGlobalMatrix(coefFunc, param, FES2T31::elementMass, tp, A, label);
    }

    public void assembleMass(double[] constCoef, BVPType tp, SRMatrix A) {
        assembleGlobalMatrix(constCoef, FES2T31::elementMass, tp, A);
    }

    public void assembleMass(double[] constCoef, BVPType tp, SCMatrix A) {
        assembleGlobalMatrix(constCoef, FES2T31::elementMass, tp, A);
    }

    public void assembleMass(double[] constCoef, BVPType tp, SRMatrix A, int... label) {
        assembleGlobalMatrix(constCoef, FES2T31::elementMass, tp, A, label);
    }

    public void assembleMass(double[] constCoef, BVPType tp, SCMatrix A, int... label) {
        assembleGlobalMatrix(constCoef, FES2T31::elementMass, tp, A, label);
    }

    public void assembleMass(double[][] vecCoef, BVPType tp, SRMatrix A) {
        assembleGlobalMatrix(vecCoef, FES2T31::elementMass, tp, A);
    }

    public void assembleMass(ScalarCFunc[] coefCFunc, Complex[] param, BVPType tp, SCMatrix A) {
        assembleGlobalMatrix(coefCFunc, param, FES2T31::elementCMass, tp, A);
    }

    public void assembleMass(Complex[] constCCoef, BVPType tp, SCMatrix A) {
        assembleGlobalMatrix(constCCoef, FES2T31::elementCMass, tp, A);
    }

    public void assembleMass(Complex[][] vecCCoef, BVPType tp, SCMatrix A) {
        assembleGlobalMatrix(vecCCoef, FES2T31::elementCMass, tp, A);
    }


    public static void elementMass(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] elemMatrix) {
        double area = CG.area2DT3(xy);
 //       System.out.println(MVO.toString(xy));
 //       System.out.println(area);
        double[] ccoef = coef[0].clone();
        MVO.fill(elemMatrix, 0.0);
        double mcoef = 0.0;
        if (tp == BVPType.CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i] *= area / 60.0;
            }
            mcoef = MVO.sum(ccoef);
            elemMatrix[0][0] = 2 * (2 * ccoef[0] + mcoef);
            elemMatrix[1][1] = 2 * (2 * ccoef[1] + mcoef);
            elemMatrix[2][2] = 2 * (2 * ccoef[2] + mcoef);
            elemMatrix[0][1] = elemMatrix[1][0] = (ccoef[0] + ccoef[1] + mcoef);
            elemMatrix[0][2] = elemMatrix[2][0] = (ccoef[0] + ccoef[2] + mcoef);
            elemMatrix[1][2] = elemMatrix[2][1] = (ccoef[1] + ccoef[2] + mcoef);

            // System.out.println(MVO.toString(elemMatrix));
        } else if (tp == BVPType.LUMP_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i] *= area / 12.0;
            }
            mcoef = MVO.sum(ccoef);
            elemMatrix[0][0] = ccoef[0] + mcoef;
            elemMatrix[1][1] = ccoef[1] + mcoef;
            elemMatrix[2][2] = ccoef[2] + mcoef;
        } else if (tp == BVPType.AXIS_CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i] *= xy[i][0] * area / 60.0;
            }
            mcoef = MVO.sum(ccoef);
            elemMatrix[0][0] = 2 * (2 * ccoef[0] + mcoef);
            elemMatrix[1][1] = 2 * (2 * ccoef[1] + mcoef);
            elemMatrix[2][2] = 2 * (2 * ccoef[2] + mcoef);
            elemMatrix[0][1] = elemMatrix[1][0] = (ccoef[0] + ccoef[1] + mcoef);
            elemMatrix[0][2] = elemMatrix[2][0] = (ccoef[0] + ccoef[2] + mcoef);
            elemMatrix[1][2] = elemMatrix[2][1] = (ccoef[1] + ccoef[2] + mcoef);
        } else if (tp == BVPType.AXIS_LUMP_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i] *= xy[i][0] * area / 12.0;
            }
            mcoef = MVO.sum(ccoef);
            elemMatrix[0][0] = ccoef[0] + mcoef;
            elemMatrix[1][1] = ccoef[1] + mcoef;
            elemMatrix[2][2] = ccoef[2] + mcoef;
        } else {
            throw new IllegalArgumentException("Wrong BVPType!");
        }
    }

    public static void elementPMLMass(int[] ele, Complex[][] xy, double[][] coef, BVPType tp, Complex[][] elemMatrix) {
        Complex area = CG.area2DT3(xy);
        double[] ccoef = coef[0].clone();
        double mcoef = 0.0;
        if (tp == BVPType.CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i] *= 1 / 60.0;
            }
            mcoef = MVO.sum(ccoef);
            elemMatrix[0][0] = area.mul(2 * (2 * ccoef[0] + mcoef));
            elemMatrix[1][1] = area.mul(2 * (2 * ccoef[1] + mcoef));
            elemMatrix[2][2] = area.mul(2 * (2 * ccoef[2] + mcoef));
            elemMatrix[0][1] = elemMatrix[1][0] = area.mul(ccoef[0] + ccoef[1] + mcoef);
            elemMatrix[0][2] = elemMatrix[2][0] = area.mul(ccoef[0] + ccoef[2] + mcoef);
            elemMatrix[1][2] = elemMatrix[2][1] = area.mul(ccoef[1] + ccoef[2] + mcoef);

            // System.out.println(MVO.toString(elemMatrix));
        }
    }

    public static void elementPMLMass(Complex[][] xy, Complex[][] coef, BVPType tp, Complex[][] elemMatrix) {
        Complex area = CG.area2DT3(xy);
        int n = coef[0].length;
        Complex[] ccoef = Complex.zeros(n);
        for (int i = 0; i < n; i++) {
            ccoef[i].real = coef[0][i].real;
            ccoef[i].image = coef[0][i].image;
        }
        Complex mcoef = new Complex(0.0, 0.0);
        if (tp == BVPType.CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i].real *= 1 / 60.0;
                ccoef[i].real *= 1 / 60.0;
            }
            mcoef = Complex.sum(ccoef);

            elemMatrix[0][0] = area.mul((ccoef[0].mul(2).add(mcoef)).mul(2));
            elemMatrix[1][1] = area.mul((ccoef[0].mul(2).add(mcoef)).mul(2));
            elemMatrix[2][2] = area.mul((ccoef[1].mul(2).add(mcoef)).mul(2));
            elemMatrix[0][1] = elemMatrix[1][0] = area.mul(ccoef[0].add(ccoef[1]).add(mcoef));
            elemMatrix[0][2] = elemMatrix[2][0] = area.mul(ccoef[0].add(ccoef[2]).add(mcoef));
            elemMatrix[1][2] = elemMatrix[2][1] = area.mul(ccoef[1].add(ccoef[2]).add(mcoef));

            // System.out.println(MVO.toString(elemMatrix));
        }
    }

    public static void elementCMass(int[] ele, double[][] xy, Complex[][] coef, BVPType tp, Complex[][] elemMatrix) {
        double area = CG.area2DT3(xy);
        Complex[] ccoef = new Complex[3];
        for (int i = 0; i < 3; i++) {
            ccoef[i] = new Complex(coef[0][i]);
        }

        Complex mcoef = null;
        if (tp == BVPType.CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i].real *= area / 60.0;
                ccoef[i].image *= area / 60.0;
            }
            mcoef = Complex.sum(ccoef);

            elemMatrix[0][0].real = 2 * (2 * ccoef[0].real + mcoef.real);
            elemMatrix[0][0].image = 2 * (2 * ccoef[0].image + mcoef.image);
            elemMatrix[1][1].real = 2 * (2 * ccoef[1].real + mcoef.real);
            elemMatrix[1][1].image = 2 * (2 * ccoef[1].image + mcoef.image);
            elemMatrix[2][2].real = 2 * (2 * ccoef[2].real + mcoef.real);
            elemMatrix[2][2].image = 2 * (2 * ccoef[2].image + mcoef.image);
            elemMatrix[0][1].real = elemMatrix[1][0].real = (ccoef[0].real + ccoef[1].real + mcoef.real);
            elemMatrix[0][1].image = elemMatrix[1][0].image = (ccoef[0].image + ccoef[1].image + mcoef.image);
            elemMatrix[0][2].real = elemMatrix[2][0].real = (ccoef[0].real + ccoef[2].real + mcoef.real);
            elemMatrix[0][2].image = elemMatrix[2][0].image = (ccoef[0].image + ccoef[2].image + mcoef.image);
            elemMatrix[1][2].real = elemMatrix[2][1].real = (ccoef[1].real + ccoef[2].real + mcoef.real);
            elemMatrix[1][2].image = elemMatrix[2][1].image = (ccoef[1].image + ccoef[2].image + mcoef.image);

//			mcoef = Complex.sum(ccoef);
//			elemMatrix[0][0].real = area / 6.0 * mcoef.real /3.0;
//			elemMatrix[0][0].image = area / 6.0 * mcoef.image /3.0;
//			elemMatrix[1][1].real = area / 6.0 * mcoef.real /3.0;
//			elemMatrix[1][1].image = area / 6.0 * mcoef.image /3.0;
//			elemMatrix[2][2].real = area / 6.0 * mcoef.real /3.0;
//			elemMatrix[2][2].image = area / 6.0 * mcoef.image /3.0;
//			elemMatrix[0][1].real = elemMatrix[1][0].real = elemMatrix[0][2].real = elemMatrix[2][0].real = elemMatrix[1][2].real = elemMatrix[2][1].real = 0.5 * elemMatrix[0][0].real;
//			elemMatrix[0][1].image = elemMatrix[1][0].image = elemMatrix[0][2].image = elemMatrix[2][0].image = elemMatrix[1][2].image = elemMatrix[2][1].image = 0.5 * elemMatrix[0][0].image;
            // System.out.println(MVO.toString(elemMatrix));
        } else if (tp == BVPType.LUMP_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i].real *= area / 12.0;
                ccoef[i].image *= area / 12.0;
            }
            mcoef = Complex.sum(ccoef);
            elemMatrix[0][0].real = ccoef[0].real + mcoef.real;
            elemMatrix[0][0].image = ccoef[0].image + mcoef.image;
            elemMatrix[1][1].real = ccoef[1].real + mcoef.real;
            elemMatrix[1][1].image = ccoef[1].image + mcoef.image;
            elemMatrix[2][2].real = ccoef[2].real + mcoef.real;
            elemMatrix[2][2].image = ccoef[2].image + mcoef.image;
        } else if (tp == BVPType.AXIS_CONSIST_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i].real *= xy[i][0] * area / 60.0;
                ccoef[i].image *= xy[i][0] * area / 60.0;
            }
            mcoef = Complex.sum(ccoef);
            elemMatrix[0][0].real = 2 * (2 * ccoef[0].real + mcoef.real);
            elemMatrix[0][0].image = 2 * (2 * ccoef[0].image + mcoef.image);
            elemMatrix[1][1].real = 2 * (2 * ccoef[1].real + mcoef.real);
            elemMatrix[1][1].image = 2 * (2 * ccoef[1].image + mcoef.image);
            elemMatrix[2][2].real = 2 * (2 * ccoef[2].real + mcoef.real);
            elemMatrix[2][2].image = 2 * (2 * ccoef[2].image + mcoef.image);
            elemMatrix[0][1].real = elemMatrix[1][0].real = (ccoef[0].real + ccoef[1].real + mcoef.real);
            elemMatrix[0][1].image = elemMatrix[1][0].image = (ccoef[0].image + ccoef[1].image + mcoef.image);
            elemMatrix[0][2].real = elemMatrix[2][0].real = (ccoef[0].real + ccoef[2].real + mcoef.real);
            elemMatrix[0][2].image = elemMatrix[2][0].image = (ccoef[0].image + ccoef[2].image + mcoef.image);
            elemMatrix[1][2].real = elemMatrix[2][1].real = (ccoef[1].real + ccoef[2].real + mcoef.real);
            elemMatrix[1][2].image = elemMatrix[2][1].image = (ccoef[1].image + ccoef[2].image + mcoef.image);
        } else if (tp == BVPType.AXIS_LUMP_MASS) {
            for (int i = 0; i < 3; i++) {
                ccoef[i].real *= xy[i][0] * area / 12.0;
                ccoef[i].image *= xy[i][0] * area / 12.0;
            }
            mcoef = Complex.sum(ccoef);
            elemMatrix[0][0].real = ccoef[0].real + mcoef.real;
            elemMatrix[0][0].image = ccoef[0].image + mcoef.image;
            elemMatrix[1][1].real = ccoef[1].real + mcoef.real;
            elemMatrix[1][1].image = ccoef[1].image + mcoef.image;
            elemMatrix[2][2].real = ccoef[2].real + mcoef.real;
            elemMatrix[2][2].image = ccoef[2].image + mcoef.image;
        } else {
            throw new IllegalArgumentException("Wrong BVPType!");
        }

    }

//	public void assembleHeatThetaScheme(ScalarCoefFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
//		double[][] kxx = new double[1][nPerEle];
//		double[][] rho = new double[1][nPerEle];
//		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
//		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
//		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
//		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
//		double dt = param[0];
//		double theta = param[1];
//		for (int i = 0; i < nE; i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coordRegion);
//			int label = mesh.getElementLabel(i);
//			for (int j = 0, nj = nPerEle; j < nj; j++) {
//				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
//				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
//			}
//			elementStiff(coordRegion, kxx, tp[0], eleMatrixStiff);
//			elementMass(coordRegion, rho, tp[1], eleMatrixMass);
//			for(int ii = 0; ii < dofPerElement; ii++ ) {
//				for(int jj = 0; jj < dofPerElement; jj++) {
//					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
//					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
//				}
//			}
//			LMatrix.assemble(ele, eleLMatrix);
//			RMatrix.assemble(ele, eleRMatrix);
//		}
//	}
//	
//	public void assembleHeatThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
//		double[][] kxx = new double[1][nPerEle];
//		double[][] rho = new double[1][nPerEle];
//		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
//		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
//		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
//		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
//		Arrays.fill(kxx[0], coefFunc[0]);
//		Arrays.fill(rho[0], coefFunc[1]);
//		double dt = coefFunc[2];
//		double theta = coefFunc[3];		
//		for (int i = 0; i < nE; i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coordRegion);
//			elementStiff(coordRegion, kxx, tp[0], eleMatrixStiff);
//			elementMass(coordRegion, rho, tp[1], eleMatrixMass);
////			System.out.println(MVO.toString(eleMatrixMass));
////			System.out.println(MVO.toString(eleMatrixStiff));
////			System.exit(0);
//			for(int ii = 0; ii < dofPerElement; ii++ ) {
//				for(int jj = 0; jj < dofPerElement; jj++) {
//					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
//					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
//				}
//			}
//			LMatrix.assemble(ele, eleLMatrix);
//			RMatrix.assemble(ele, eleRMatrix);
//		}
//	}

    public void assembleThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix,
                                    SRMatrix RMatrix) {
        assembleHeatThetaScheme0(coefFunc, param, tp, FES2T31::elementStiff, FES2T31::elementMass, LMatrix, RMatrix);
    }

    public void assembleThetaScheme(double[] coefFunc, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
        assembleHeatThetaScheme0(coefFunc, tp, FES2T31::elementStiff, FES2T31::elementMass, LMatrix, RMatrix);
    }



    public static void elementSource(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
        double area = CG.area2DT3(xy);
        double sumF = MVO.sum(f[0]);
        Arrays.fill(eleF, 0.0);
        switch (tp) {
            case COMMON:
            case CONSTITUTE_MATRIX_COMMON:
                for (int i = 0; i < 3; i++) {
                    eleF[i] = (f[0][i] + sumF) * area / 12.0;
                }
                break;
            case CONSTITUTE_MATRIX_AXIS:
            case AXISYMMETRIC:
                for (int j = 0; j < 4; j++) {
                    for (int i = 0; i < 3; i++) {
                        eleF[i] += area * TW3[j] * MVO.dot_product(TP3[j], f[0])
                                * MVO.dot_product(TP3[j], new double[]{xy[0][0], xy[1][0], xy[2][0]}) * TP3[j][i];
                    }
                }
                break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
    }

    public void assembleConvection(double[] constCoef, BVPType tp, SRMatrix A) {
        assembleGlobalMatrix(constCoef, FES2T31::elementConvection, tp, A);
    }

    public static void elementConvection(int[] ele, double[][] xy, double[][] v, BVPType tp, double[][] eleMatrix) {

        MVO.fill(eleMatrix, 0.0);
        double[][] grad = new double[3][2];
        double area = CG.gradient2DT3(xy, grad);
        area /= 12.0;

        double v11 = v[0][0];
        double v12 = v[0][1];
        double v13 = v[0][2];

        double v21 = v[1][0];
        double v22 = v[1][1];
        double v23 = v[1][2];

        double v1N1 = (v11 * 2 + v12 + v13);
        double v2N1 = (v21 * 2 + v22 + v23);

        double v1N2 = (v11 + v12 * 2 + v13);
        double v2N2 = (v21 + v22 * 2 + v23);

        double v1N3 = (v11 + v12 + v13 * 2);
        double v2N3 = (v21 + v22 + v23 * 2);

        for (int i = 0; i < 3; i++) {
            eleMatrix[0][i] = (v1N1 * grad[i][0] + v2N1 * grad[i][1]) * area;
            eleMatrix[1][i] = (v1N2 * grad[i][0] + v2N2 * grad[i][1]) * area;
            eleMatrix[2][i] = (v1N3 * grad[i][0] + v2N3 * grad[i][1]) * area;
        }
        // System.out.println(MVO.toString(eleMatrix));
    }

    public static void elementFlux(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] k) {
        double len = CG.length2DL2(coord);
        double sumG = MVO.sum(g[0]);
        double sumX = coord[0][0] + coord[1][0];
        switch (tp) {
            case COMMON:
                for (int i = 0; i < 2; i++) {
                    k[i] = (sumG + g[0][i]) * len / 6.0;
                }
                break;
            case AXISYMMETRIC:
                for (int i = 0; i < 2; i++) {
                    k[i] = (sumG * sumX + 2 * g[0][i] * coord[i][0]) * len / 12.0;
                }
                break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
    }

    public static void elementCFlux(int[] ele, double[][] coord, Complex[][] g, BVPType tp, Complex[] k) {
        double len = CG.length2DL2(coord);
        Complex sum = MVO.sum(g[0]);
        double sumX = coord[0][0] + coord[1][0];
        switch (tp) {
            case COMMON:
                for (int i = 0; i < 2; i++) {
                    k[i].real = (sum.real + g[0][i].real) * len / 6.0;
                    k[i].image = (sum.image + g[0][i].image) * len / 6.0;
                }
                break;
            case AXISYMMETRIC:
                for (int i = 0; i < 2; i++) {
                    k[i].real = (sum.real * sumX + 2 * g[0][i].real * coord[i][0]) * len / 12.0;
                    k[i].image = (sum.image * sumX + 2 * g[0][i].image * coord[i][0]) * len / 12.0;
                }
                break;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
    }


    public void assembleBoundaryMass(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(coefFunc, param, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void assembleBoundaryMass(double[] constConvect, BVPType tp, SRMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(constConvect, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void assembleBoundaryMass(double[][] vectorConvect, BVPType tp, SRMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(vectorConvect, FES2T31::elementBoundaryMass, tp, A, bdLabel);
    }

    public void assembleBoundaryMass(ScalarCFunc[] coefFunc, Complex[] param, BVPType tp, SCMatrix A, int... bdLabel) {
        assembleBoundaryMatrix(coefFunc, param, FES2T31::elementCBoundaryMass, tp, A, bdLabel);
    }

    public static void elementBoundaryMass(int[] ele, double[][] xy, double[][] h, BVPType tp, double[][] eleMatrix) {
        double len = CG.length2DL2(xy);
        double[] x = MVO.col(xy, 0);
        double[] g = h[0];
        MVO.fill(eleMatrix, 0.0);
        switch (tp) {
            case CONSIST_MASS:
                len /= 12;
                eleMatrix[0][0] = (3 * g[0] + g[1]) * len;
                eleMatrix[0][1] = (g[0] + g[1]) * len;
                eleMatrix[1][0] = eleMatrix[0][1];
                eleMatrix[1][1] = (g[0] + 3 * g[1]) * len;
                break;
            case LUMP_MASS:
                len /= 6;
                eleMatrix[0][0] = (2 * g[0] + g[1]) * len;
                eleMatrix[1][1] = (g[0] + 2 * g[1]) * len;
                break;
            case AXIS_CONSIST_MASS:
                len /= 60;
                eleMatrix[0][0] = (12 * g[0] * x[0] + 3 * (g[1] * x[0] + g[0] * x[1]) + 2 * g[1] * x[1]) * len;
                eleMatrix[0][1] = (3 * (g[0] * x[0] + g[1] * x[1]) + 2 * (g[0] * x[1] + g[1] * x[0])) * len;
                eleMatrix[1][0] = eleMatrix[0][1];
                eleMatrix[1][1] = (2 * g[0] * x[0] + 3 * (g[1] * x[0] + g[0] * x[2]) + 12 * g[1] * x[1]) * len;
            case AXIS_LUMP_MASS:
                len /= 12;
                eleMatrix[0][0] = (2 * g[0] * x[0] + MVO.sum(g) * MVO.sum(x)) * len;
                eleMatrix[1][1] = (2 * g[1] * x[1] + MVO.sum(g) * MVO.sum(x)) * len;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
    }

    public static void elementCBoundaryMass(int[] ele, double[][] xy, Complex[][] h, BVPType tp, Complex[][] eleMatrix) {
        double len = CG.length2DL2(xy);
        double[] x = MVO.col(xy, 0);
        Complex[] g = h[0];
        switch (tp) {
            case CONSIST_MASS:
                len /= 12;
                eleMatrix[0][0].real = (3 * g[0].real + g[1].real) * len;
                eleMatrix[0][0].image = (3 * g[0].image + g[1].image) * len;
                eleMatrix[0][1].real = (g[0].real + g[1].real) * len;
                eleMatrix[0][1].image = (g[0].image + g[1].image) * len;
                eleMatrix[1][0].real = eleMatrix[0][1].real;
                eleMatrix[1][0].image = eleMatrix[0][1].image;
                eleMatrix[1][1].real = (g[0].real + 3 * g[1].real) * len;
                eleMatrix[1][1].image = (g[0].image + 3 * g[1].image) * len;
                break;
            case LUMP_MASS:
                len /= 6;
                eleMatrix[0][0].real = (2 * g[0].real + g[1].real) * len;
                eleMatrix[1][1].real = (g[0].real + 2 * g[1].real) * len;
                eleMatrix[0][0].image = (2 * g[0].image + g[1].image) * len;
                eleMatrix[1][1].image = (g[0].image + 2 * g[1].image) * len;
                break;
            case AXIS_CONSIST_MASS:
                len /= 60;
                eleMatrix[0][0].real = (12 * g[0].real * x[0] + 3 * (g[1].real * x[0] + g[0].real * x[1]) + 2 * g[1].real * x[1]) * len;
                eleMatrix[0][1].real = (3 * (g[0].real * x[0] + g[1].real * x[1]) + 2 * (g[0].real * x[1] + g[1].real * x[0])) * len;
                eleMatrix[1][0].real = eleMatrix[0][1].real;
                eleMatrix[1][1].real = (2 * g[0].real * x[0] + 3 * (g[1].real * x[0] + g[0].real * x[2]) + 12 * g[1].real * x[1]) * len;
                eleMatrix[0][0].image = (12 * g[0].image * x[0] + 3 * (g[1].image * x[0] + g[0].image * x[1]) + 2 * g[1].image * x[1]) * len;
                eleMatrix[0][1].image = (3 * (g[0].image * x[0] + g[1].image * x[1]) + 2 * (g[0].image * x[1] + g[1].image * x[0])) * len;
                eleMatrix[1][0].image = eleMatrix[0][1].image;
                eleMatrix[1][1].image = (2 * g[0].image * x[0] + 3 * (g[1].image * x[0] + g[0].image * x[2]) + 12 * g[1].image * x[1]) * len;
            case AXIS_LUMP_MASS:
                len /= 12;
                eleMatrix[0][0].real = (2 * g[0].real * x[0] + MVO.sum(g).real * MVO.sum(x)) * len;
                eleMatrix[1][1].real = (2 * g[1].real * x[1] + MVO.sum(g).real * MVO.sum(x)) * len;
                eleMatrix[0][0].image = (2 * g[0].image * x[0] + MVO.sum(g).image * MVO.sum(x)) * len;
                eleMatrix[1][1].image = (2 * g[1].image * x[1] + MVO.sum(g).image * MVO.sum(x)) * len;
            default:
                throw new IllegalArgumentException("Wrong BVPType!");
        }
    }

    public double computeL2Norm(double[] u) {
        double result = 0.0;
        for (int i = 0; i < mesh.getNt(); i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            double[] uu = MVO.getSlice(u, ele);
            result += area
                    * (uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2] + uu[0] * uu[1] + uu[1] * uu[2] + uu[2] * uu[0])
                    / 6;
        }
        return Math.sqrt(result);
    }

    public double computeH1SemiNorm(double[] u) {
        double result = 0.0;
        double[][] grad = new double[3][2];
        double dxU, dyU;
        for (int i = 0; i < mesh.getNt(); i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            CG.gradient2DT3(coordRegion, grad);
            double[] uu = MVO.getSlice(u, ele);
            double[] dx = MVO.col(grad, 0);
            double[] dy = MVO.col(grad, 1);
            dxU = MVO.dot_product(dx, uu);
            dyU = MVO.dot_product(dy, uu);
            result = result + area * (dxU * dxU + dyU * dyU);
        }
        return Math.sqrt(result);
    }

    public void computeGradient(double[] u, double[][] gradU) {
        int len = u.length;
        ParamCheck.checkEqual(len, this.ndof);
        ParamCheck.checkEqual(len, gradU[0].length);
        ParamCheck.checkEqual(len, gradU[1].length);

        double[] angle = new double[3];
        double[][] coord = new double[3][2];
        double[][] grad = new double[3][2];
        double[] angleCount = new double[this.GetNdof()];
        double dxValue, dyValue;
        double[] eleU = new double[3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coord);
            CG.gradient2DT3(coord, grad);
            CG.angle2DT3(coord, angle);
            eleU[0] = u[ele[0]];
            eleU[1] = u[ele[1]];
            eleU[2] = u[ele[2]];
            double[] dx = MVO.col(grad, 0);
            double[] dy = MVO.col(grad, 1);
            dxValue = MVO.dot_product(eleU, dx);
            dyValue = MVO.dot_product(eleU, dy);
            for (int j = 0; j < 3; j++) {
                gradU[0][ele[j]] += dxValue * angle[j];
                gradU[1][ele[j]] += dyValue * angle[j];
                angleCount[ele[j]] += angle[j];
            }
        }
        for (int i = 0, n = mesh.getNv(); i < n; i++) {
            gradU[0][i] /= angleCount[i];
            gradU[1][i] /= angleCount[i];
        }
    }

    /**
     *
     * @param u
     * @param gradU
     */
    public void computeGradientBySPR(double[] u, double[][] gradU) {
        int len = u.length;
        ParamCheck.checkEqual(len, this.ndof);
        ParamCheck.checkEqual(len, gradU[0].length);
        ParamCheck.checkEqual(len, gradU[1].length);

        double[][] coord = new double[3][2];
        double[][] grad = new double[3][2];
        double[] eleU = new double[3];
        double[][] elementGradient = new double[2][mesh.getNt()];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coord);
            CG.gradient2DT3(coord, grad);
            eleU[0] = u[ele[0]];
            eleU[1] = u[ele[1]];
            eleU[2] = u[ele[2]];
            double[] dx = MVO.col(grad, 0);
            double[] dy = MVO.col(grad, 1);
            elementGradient[0][i] = MVO.dot_product(eleU, dx);
            elementGradient[1][i] = MVO.dot_product(eleU, dy);
        }
        System.out.println("Element Gradient: ");
        System.out.println("dx min: " + MVO.min(elementGradient[0]) +" max: " + MVO.max(elementGradient[0]));
        System.out.println("dy min: " + MVO.min(elementGradient[1]) +" max: " + MVO.max(elementGradient[1]));
        HashSet<Integer>[] patchSet = mesh.SuperconvergencePatches(3);
        int nvIdx = 0;
        for(HashSet<Integer> patch: patchSet) {
            System.out.println(patch);
            int pSize = patch.size();
            double[][] A = new double[pSize][];
            double[] bx = new double[pSize];
            double[] by = new double[pSize];
            int pIdx = 0;
            for(int idx : patch) {
                int[] ele = mesh.getElement(idx);
                mesh.getCoordPerElement(ele, coord);
                double[] xy = mesh.getinteriorPointInElement(coord);
                A[pIdx] = new double[]{1.0, xy[0], xy[1]};
                bx[pIdx] = elementGradient[0][idx];
                by[pIdx] = elementGradient[1][idx];
                pIdx++;
            }
            double[][] AT = MVO.transpose(A);
            double[][] ATA = MVO.matmul(AT, A);
            double[] dx = new double[3];
            double[] dy = new double[3];
            double[] ATBx = MVO.matmul(AT, bx);
            double[] ATBy = MVO.matmul(AT, by);
            DirectDSolver solver = new DirectDSolver(ATA);
            solver.CholeskyFactorize();
            solver.CholeskySolve(ATBx, dx);
            solver.CholeskySolve(ATBy, dy);
            double[] xy = mesh.getNode(nvIdx);
            double[] P = new double[]{1.0, xy[0], xy[1]};
            gradU[0][nvIdx] = MVO.dot_product(P, dx);
            gradU[1][nvIdx] = MVO.dot_product(P, dy);
            nvIdx++;
        }
    }

    public double[] computeErrorIndictor_MATLAB(double[] u, ScalarFunc[] coefFunc, double[] params) {
        double[] eleError = new double[mesh.getNt()];
        Mesh2T3 mesh = (Mesh2T3) (this.mesh);
        SIMatrix edgeMat = mesh.edgeFromElement();
        double[][] coef = new double[coefFunc.length][nPerEle];
        double alpha = params[0];
        double beta = params[1];
        double m = params[2];
        double[][] grad = new double[3][2];
        double[] L = new double[3];
        double[] nx = new double[3];
        double[] ny = new double[3];
        double[] eleU = new double[3];
        double[] err1 = new double[3];
        double[] dx = null;
        double[] dy = null;
        double dxU, dyU, ndxU, ndyU;
        double maxL = 0;
        double area;
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            int label = mesh.getElementLabel(i);
            getCoefFromFunc(coefFunc, coordRegion, label, null, coef);
            double[] kxx = coef[0];
            double[] aa = coef[1];
            double[] source = coef[2];
            area = CG.gradient2DT3(coordRegion, grad);
            CG.edge2DT3(coordRegion, L, nx, ny);
            for (int j = 0; j < 3; j++) {
                eleU[j] = u[ele[j]];
            }
            maxL = MVO.max(L);
            dx = MVO.col(grad, 0);
            dy = MVO.col(grad, 1);
            dxU = MVO.dot_product(eleU, dx);
            dyU = MVO.dot_product(eleU, dy);
            double mkxx = MVO.average(kxx);
            for (int j = 0; j < 3; j++) {
                err1[j] = source[j] - aa[j] * eleU[j];
            }
            double temp = area * (err1[0] * err1[0] + err1[1] * err1[1] + err1[2] * err1[2] + err1[0] * err1[1]
                    + err1[1] * err1[2] + err1[2] * err1[0]) / 6.0;
            double errorOnArea = alpha * Math.pow(maxL, m) * Math.sqrt(temp);

            double jump1 = 0.0;
            int neighEle = edgeMat.getElement(ele[1], ele[0]);
            if (neighEle != 0) {
                int[] neighElement = mesh.getElement(neighEle);
                mesh.getCoordPerElement(neighElement, coordRegion);
                CG.gradient2DT3(coordRegion, grad);
                dx = MVO.col(grad, 0);
                dy = MVO.col(grad, 1);
                for (int j = 0; j < 3; j++) {
                    eleU[j] = u[neighElement[j]];
                }
                ndxU = MVO.dot_product(eleU, dx);
                ndyU = MVO.dot_product(eleU, dy);
                jump1 = ((dxU - ndxU) * nx[0] + (dyU - ndyU) * dy[0]) * mkxx;
                jump1 = 0.50 * Math.pow(L[0], 2 * m) * jump1 * jump1;
            }

            double jump2 = 0.0;
            neighEle = edgeMat.getElement(ele[2], ele[1]);
            if (neighEle != 0) {
                int[] neighElement = mesh.getElement(neighEle);
                mesh.getCoordPerElement(neighElement, coordRegion);
                CG.gradient2DT3(coordRegion, grad);
                dx = MVO.col(grad, 0);
                dy = MVO.col(grad, 1);
                for (int j = 0; j < 3; j++) {
                    eleU[j] = u[neighElement[j]];
                }
                ndxU = MVO.dot_product(eleU, dx);
                ndyU = MVO.dot_product(eleU, dy);
                jump2 = ((dxU - ndxU) * nx[1] + (dyU - ndyU) * dy[1]) * mkxx;
                jump2 = 0.50 * Math.pow(L[1], 2 * m) * jump2 * jump2;
            }

            double jump3 = 0.0;
            neighEle = edgeMat.getElement(ele[0], ele[2]);
            if (neighEle != 0) {
                int[] neighElement = mesh.getElement(neighEle);
                mesh.getCoordPerElement(neighElement, coordRegion);
                CG.gradient2DT3(coordRegion, grad);
                dx = MVO.col(grad, 0);
                dy = MVO.col(grad, 1);
                for (int j = 0; j < 3; j++) {
                    eleU[j] = u[neighElement[j]];
                }
                ndxU = MVO.dot_product(eleU, dx);
                ndyU = MVO.dot_product(eleU, dy);
                jump3 = ((dxU - ndxU) * nx[2] + (dyU - ndyU) * dy[2]) * mkxx;
                jump3 = 0.50 * Math.pow(L[2], 2 * m) * jump3 * jump3;
            }
            double errorOnLine = beta * Math.sqrt(jump1 + jump2 + jump3);
            eleError[i] = errorOnArea + errorOnLine;
        }
        return eleError;
    }

    public double[] integrate(ScalarFunc[] sfun, double[] param) {
        int nf = sfun.length;
        double[] result = new double[nf];
        double[][] coef = new double[nf][3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            int label = mesh.getElementLabel(i);
            double area = CG.area2DT3(coordRegion);
            getCoefFromFunc(sfun, coordRegion, label, param, coef);
            for (int j = 0; j < nf; j++) {
                result[j] += MVO.sum(coef[j]) * area / 3.0;
            }

        }
        return result;
    }

    public double integrate(double[] sol) {
        double result = 0.0;
        double[] coef = new double[3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            coef[0] = sol[ele[0]];
            coef[1] = sol[ele[1]];
            coef[2] = sol[ele[2]];
            result += MVO.sum(coef) * area / 3.0;
        }
        return result;
    }

    public double[] integrate(double[][] sol) {
        int nf = sol.length;
        double[] result = new double[nf];
        double[][] coef = new double[nf][3];
        for (int i = 0, n = mesh.getNt(); i < n; i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            for (int j = 0; j < nf; j++) {
                coef[j][0] = sol[j][ele[0]];
                coef[j][1] = sol[j][ele[1]];
                coef[j][2] = sol[j][ele[2]];
                result[j] += MVO.sum(coef[j]) * area / 3.0;
            }
        }

        return result;
    }

    public void map(double[] source, FES fs, double[] destination) {
        int nv1 = mesh.getNv();
        int nv2 = fs.getMesh().getNv();
        Mesh mesh2 = fs.getMesh();
        ParamCheck.checkEqual(nv1, source.length);
        ParamCheck.checkEqual(nv2, destination.length);
        for (int i = 0; i < nv2; i++) {
            double[] xy = mesh2.getNode(i);
            double interp = getValue(source, xy);
            if (interp != 1.0e30) {
                destination[i] = interp;
            }
        }
    }

}
